Stable multispeed lattice Boltzmann methods 
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We demonstrate how to produce a stable multispeed lattice Boltzmann method (LBM) for a wide 
range of velocity sets, many of which were previously thought to be intrinsically unstable. We use 
non-Gauss-Hermitian cubatures. The method operates stably for almost zero viscosity, has second- 
order accuracy, suppresses typical spurious oscillation (only a modest Gibbs effect is present) and 
introduces no artificial viscosity. There is almost no computational cost for this innovation. 

DISCLAIMER: Additional tests and wide discussion of this preprint show that the 
claimed property of coupled steps: no artificial dissipation and the second-order ac- 
curacy of the method are valid only on sufficiently fine grids. For coarse grids the 
higher-order terms destroy coupling of steps and additional dissipation appears. 

The equations are true. 



I. INTRODUCTION 

The lattice Boltzmann method (LBM) is a discrete ve- 
locity method primarily used in the numerical simulation 
of complex fluid problems. The essence of the method is 
that a finite number of populations stream and collide on 
a fixed computational lattice in such a manner that their 
populations' velocity moments obey the Navier-Stokes 
equations. 

There are a number of ways to derive the method 
Historically, the method comes from lattice gas automata 
theory but it can also be derived directly by the overre- 
laxation discretization of Boltzmann's kinetic transport 
equation. 

We prefer to describe the LBM without direct reference 
to Boltzmann's equation. Indeed, we prefer to describe 
the LBM as being generated wholly by the mechanical 
motion, of a single-particle distribution function, by free- 
flight and entropic involution. In this setting, the LBM is 
the discrete dynamical system which arises by discretiz- 
ing this motion in a particular manner. The discrete 
velocity set arises as approximation nodes of certain cu- 
bature in velocity space 0, where these nodes are auto- 
morphisms of some underlying lattice. 

A consequence of this realisation is that the stability 
analysis (and analysis of conservation laws) of the LBM 
is more natural. One is able to clearly identify the main 
instability mechanisms of the LBM (Sect.|TT|, which are 
triggered when the continuous mechanical motion of free- 
flight and eutopic involution is discretized. 

A further consequence of viewing the LBM in this set- 
ting (also in Sect. [TTI is that a prescription of coupled 
steps is suggested to stabilise the method. After proceed- 
ing for one step of the usual overrelaxation LBM scheme 
(LBGK), populations are streamed and then replaced 
with their local equilibrium distribution. Additional dis- 
sipation results but the scheme retains the second-order 
in time accuracy of LBGK. Compared with the standard 
LBM, the proposed scheme of coupled steps constitutes 



no additional computational cost. 

Each cubature rule gives rise to its own LBM, and 
multispeed lattice schemes can be readily concocted 
(Sect.HlH). Multispeed lattices are an attractive prospect 
because they promise greater flexibility in the attain- 
able sound speed provided by the model, allow for the 
modelling of nonisothermic flows and the dynamics of 
higher moments as well. The literature sometimes con- 
tains statements whose sentiment is that some multi- 
speed lattice LBM schemes are inherently numerically 
unstable (see, e.g., On the contrary, we observe 

that the aforementioned LBM of coupled steps is sta- 
ble for a variety of multispeed lattices and we demon- 
strate this through the numerical simulation of a one- 
dimensional isothermal shock tube (Sect. ITV)). 



II. THE LATTICE BOLTZMANN METHOD 

To describe the LBM as being generated by free-flight 
and entropic involution requires a certain amount of 
background knowledge and terminology which we now 
briefly introduce. For proofs and further justification of 
any statements, the reader should consult [1, Q. The 
free-flight equation is given by 
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where / = /{x, v, t) is a single-particle distribution func- 
tion, X is the space vector, v is velocity. This equation 
preserves the value of a strictly concave entropy func- 
tional, S{f). The default choice of entropy is 



S{f) 



flogfdvdx. 



(2) 



In addition, we have a flxed linear mapping m : f t—^ AI 
to some macroscopic variables. For example, in hydrody- 
namic applications M is the vector of flve hydrodynamic 
fields (n,nu,E) (density-momentum-energy): 
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For each M , the quasiequilihrium state f^j is defined as 
the unique solution of the constrained optimisation prob- 
lem: argmax{5'(/) : m{f) = AI}. For the hydrodynamic 
fields M = M{x,t) the quasiequilihrium state is the 
well known local Maxwellian 

where m is particle mass, is Boltzmann's constant and 
T is kinetic temperature. 

For every /, a corresponding quasiequilihrium state 
f^(^fj is defined. The set of all quasiequilihrium states 
is parameterised by the macroscopic variables, M, and 
defines the quasiequilihrium manifold, which we denote 
by qo- 

The quasiequilihrium approximation for the free-flight 
equation ([T]) is the following evolution equation for the 
macroscopic variables 

^ = -m(^.V(/;,))). (5) 

For the hydrodynamics fields ([3]) the system ([5]) is the 
compressible Euler equations. 

We denote by Ilg the projector of a point / onto qg: 
: / I— > f*n(f)- Let Qt be the time shift transfor- 
mation for the free-flight equation ([1]): 8( : f{x,v) i-^ 
f{x — vt,v). Now, for a flxed time step r, the follow- 
ing step provides a second-order in time step t approx- 
imation to the solution of the conservative macroscopic 
equations ([5]): 

M(o) = m(n5(e_,/2(/M))) 

^ miUsiQrMrM))) = A/(t). (6) 

If we would like to model dissipative dynamics then we 
should follow the free-flight trajectory by some extra time 
<^ < t: 

M(o) = m(ns(e_^/2(/X,))) 

^ m(n5(e,+^/2(/Xf))) = M{, + i9)= M{t), (7) 

where = r— ^. Now, ([7]) provides a second-order in time 
step T approximation to the solution of the compressible 
Navier-Stokes equations with dynamic viscosity fi ~ 
and Prandtl number Pr = 1. 

Equations ^ and ([7]) constitute one step of an ap- 
proximation to some conservative and dissipative macro- 
scopic equations, respectively. To iterate and produce 
a step wise approximation to the macroscopic equations 
for time greater than r we can employ (partial) entropic 
involution. 

It will be useful to define the film of nonequilibrium, 
states 0, 0, [1] as the manifold q that is the trajectory 
(forward and backward in time) of the quasiequilihrium 
manifold. A point / G q is naturally parameterised by 
(M, t): / = qM,T, where M = m{f) is the value of the 
macroscopic variables, and r — T{f) is the time shift from 



a quasiequilihrium state: Q-rif) is a quasiequilihrium 
state for some (other) value of M. The quasiequilihrium 
manifold divides q into two parts, q = q_ UqoUq+, where 
q- = Ualt\ t < 0} and q+ = {gM.rl t > 0}. 

Now, the partial entropic involution operator 1^ is de- 
fined as follows: for / G q± the point g — Ig{f) G q;p, 
for (3 £ [1/2, 1), is specified by the two conditions: 

m{g) = m{f), 
S{g) - S{ns{f)) = (2/3 - 1)2(5(/) - S{ns{f))). 

If /3 = 1 in this definition we refer to the operator as 
entropic involution and denote it by Is- The point Igif), 
fi G [1/2, 1), is closer to the quasiequilihrium point n5(/) 
(with respect to entropy) than Is{f)- 

If /o S qo, then, after some initial steps, the following 
sequence gives a second-order in time step r approxi- 
mation of the compressible Navier-Stokes equations (as 
mentioned above) with = (1 — /3)t//3, f3 € [1/2, 1): 

M(nr) =TO((/fe,)"/o), 71 = 0,1,.... (8) 

For entropic involution (/3 = 1), ([5]) provides a second- 
order in time step r approximation of the compressible 
Euler equations with time step 2t. 

Finally, to generate the standard LBM we must per- 
form three tasks: 

1. transfer to a finite number of velocities with the 
same macroscopic equations; 

2. transfer from space to a lattice, where these veloc- 
ities are automorphisms; 

3. transfer from dynamics and involution on q to the 
whole space of states. 

The first two points will be addressed in Sect. Ill CI and 
Sect, mil To address the final point we replace the oper- 
ator Ig with the transformation 

/o^:/->n5(/)-l-(2/3-l)(n5(/)-/). (9) 

If, for a given /o € qo, the sequence ^ gives the sought 
after second-order in time step r approximation of the 
macroscopic equations, then the sequence 

M(nT) =TO((/o^e,)"/o), n = 0,l,..., (10) 

also gives a second-order approximation to the same 
equations. 

A. Instability mechanisms 

There are a number of sources of instability which are 
triggered when the operator Ig is replaced with /q 
Indeed, we identify three main instability mechanisms. 

Firstly, a small shift in / in the direction of the vector 
f^^sif) does not relax back for (3=1, and its relaxation 



FIG. 1: Neutral stability and one-step oscillations in a se- 
quence successive entropic involutions. Bold dotted line - a 
perturbed motion, A - direction of neutral stability. 

is very slow for /? ~ 1 (i.e., for small viscosity). This effect 
is most easily illustrated by Fig. [1] where a perturbed 
chain of entropic involutions is shown. This instability, 
which we refer to as neutral instability, causes one step 
oscillations to be triggered. 

Secondly, there is a nonlinear instability due to the 
nonlinear nature of the projector Us- More specifically, 
for all of the accuracy estimates in the previous section, 
we implicitly use the assumption that / is sufficiently 
close to qo. The sequences @ and ([TU|) give the same ac- 
curacy as one another, but a long chain of steps with the 
linearised operator Iq ^ can lead far from the quasiequi- 
librium manifold qo and even from q itself. 

Finally, there is a directional instability that can af- 
fect accuracy. There can be large deviation in the angle 
the vector / — Hsif) makes with the tangent space to 
q. The directional instability changes the structure of 
dissipation terms in the target macroscopic equations: 
accuracy is decreased to the first-order in r and signif- 
icant fluctuations of the Prandtl number and viscosity 
coefficient may occur. 

B. Stabilisation 

There is a strikingly simple prescription that simulta- 
neously alleviates the effect of all three instabilities iden- 
tified in the previous section. 

Consider starting from a point /o G qo, then evolve the 
state by Qr and apply the operator Iq (l9|), as is usual. 
Then, evolve by 8^ again and project back on to qo using 

M{0) - m(/o) ^ m(ns(e,(/o^(e,(/o))))) = M{2t). 

(11) 

This coupled step with quasiequilibrium ends is illus- 
trated in Fig. [3 and gives a second-order in time r ap- 
proximation, to the shift in time 2r, for the target macro- 
scopic equations with ^ = 2(1 — /3)r, f3 G [1/2, 1]. The 
procedure of periodically restarting from the quasiequi- 
librium manifold introduces additional dissipation of or- 
der , and the perturbation of accuracy is of order . 
Hence, the method has the second-order accuracy. 

The user should take into account that <^ significantly 
depends on the chain construction. For the sequence ([5]) 



FIG. 2: The scheme of coupled steps with quasiequilibrium 
ends (fTT|) . 



we have <r = (1 — (3)t/ (3, and for the sequence of steps 
PT|) we have <r = 2(1 — /3)t. The viscosity coefficient is 
proportional to If 1 — /3 is small, the coupled step pT|) 
gives around two times larger viscosity than ([8]). 

The result of restarting from qo is that the afore- 
mentioned neutral instability is obliterated and the ef- 
fects of nonlinear and directional instabilities are entirely 
marginalised. 



C. Entropy, energy and equilibrium 

Many specific forms of entropy for LBMs have been dis- 
cussed in the literature. There exist two methods of con- 
struction: from the Boltzmann entropy approximation to 
equilibrium (for given macroscopic variables), and from 
the equilibrium approximation to the Boltzmann entropy 
approximation. For the latter, the universal entropy for- 
mula for a discrete distribution / = {fi) is the KuUback 
entropy 

i \J7n{f),i/ 

where /^j is the quasiequilibrium distribution param- 
eterised by M. It is trivial to check that f^j — 
argmax{5'_R-(/) : m{f) = M}. The entropic involution 
conserves the values of all elements of M. Hence, for 
energy conservation it is sufficient to have energy inside 
M (the distributions / and f^(f) should have the same 
energy). In this sense, the long standing problem of en- 
ergy conservation in LBMs has an obvious and physically 
meaningful solution (see an example of such a solution 
in [l^ ) . The classical physical sense of ([12]) is a nonequi- 
librium extension of a Massieu-Planck-Kramers function 
(or grand canonical potential) for a given set of macro- 
scopic variables M. 

A constructive approach to choose discrete S and flj 
is: we use the classical (continuous) quasi-equilibrium 
distributions and cubature rules to get the appropriate 
discrete approximation. Then, the entropy that we need 
for self-consistency of our construction is given by . 
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III. CUBATURE AND MULTISPEED LATTICES 

To complete the realisation of the LBM, configuration 
space is discretized by selecting a finite number of veloc- 
ities {vi,V2, ■ ■ ■ ,vg}. More specifically, for a fixed time 
step r, the space discretization will consist of a lattice £ 
for which the velocities v^t are automorphisms. Choos- 
ing to discretize in this manner means that no spatial 
discretization error is committed when the continuous 
free-flight equations 



Method Order 



dt 



V/ = 0, z = 



are integrated for time r and evaluated on £. Here, fi 
denotes a single-particle distribution function associated 
with the velocity Vi. Note that all of the discretization 
error is contained in the selection of the discrete veloci- 
ties. 

For hydrodynamics it remains to evaluate the inte- 
gral moments ([3]). By construction, the integrals ([3]) re- 
main unchanged under the substitution of / with the 
quasiequilibrium state Our approach will be to 

construct cubature rules based upon reproduction of low 
degree polynomials using the discrete velocities, Vi, as ab- 
scissas. We replace the quasiequilibrium state (|4]) with its 
second-order Taylor expansion in u with terms involving 
and higher disregarded. We denote this polynomial 
by p^j — p*j^,f{v), and demand that its first few moments 
are evaluated exactly by the cubature rule. Another ap- 
proach, which we do not pursue here, is to employ a 
discrete version of the entropy ([2]) from the outset (see, 
e.g., i). 

Let us concern ourselves with modelling the isothermal 
Navier-Stokes equations for the remainder of the paper 
and set = 2kBT/m. The "sound speed" of the model 
is Cs where = For the isothermal model we 

should at least consider all moments up to the third- 
order, which is equivalent to finding the cubature weights 
Wi such that 

I- 2 ^ 2 

/ «'exp(-^)dt; = ^7«,i;,'=exp(-^), 

i—l 

for fc = 0, 1, 2, 3. Using this cubature, the discrete hydro- 
dynamic moments are given by the expressions 

i i i 

where fi = fi{x, t) := Wif{x, Vi,t). The LBM realisation 
of ^ is then 

Mx + v,r,t + r) - /* + (2/3 - 1)(/* - /,), (13) 

where /* = f*{x,t) :— Wip\j{vi). The name often given 
to in the literature is LBGK. The prescribed scheme 
of coupled steps pT|) is: 



fi{x+ViT,t+T) = 



TVstop odd. 



/; + (2/3-l)(/*-/,), otherwise, 

(14) 



G-H 
N-C 
N-C 
N-C 
N-C 
N-C 
N-C 
N-C 
N-C 



{0,±1} 
{0,±1,±2} 
{0,±1,±2} 
{0,±1,±3} 
{0,±1,±3} 
{0,±1,±4} 
{0,±1,±4} 



I 3' 6 J 

I 16 ' 24 ' 96-1 

I 2 ' 6 ' 12-1 
r 19 15 1 1 
t 36 ' 64 ' 576 J 
Si 3 J_l 
l 9 ' 8 ' 72-1 
r 33 _29_ 1 1 
l 64 ' 120 ' 1920 > 



1/V3 
1/72 
1 

l/%/2 
1 

1 



L 64 ' 120 ' 1920 J 
I 8 ' 30 ' 240-1 

{0,±1,±2,±3} 3I0 , l/\/2 

{0,±1,±2,±3} {^,1,^,^} 1 



TABLE L Various ID lattices and weights using difler- 
ent quadratures: G-H means Gauss-Hermite, N-C means 
Newton-Cotes. The second column shows the order of the 
moment which is exactly evaluated by the quadrature rule. 
The final column gives the lattice "sound speed". 



where Nstcp is the cumulative number of time steps taken 
in the simulation. 

Clearly, each set {vi,Wi) of velocities and associated 
cubature weights will define its own LBM. 

One way to handle the multivariate case is to use ten- 
sor and algebraic products of one-dimensional velocities 
and weights, respectively. Let us therefore restrict our 
discussion to one space dimension. 

One particular choice of quadrature is Gauss-Hermite 
quadrature. Indeed, the three-point Gauss-Hermite rule 
(which exactly evaluates the fifth order moment) is pop- 
ular [ll|, [l2|. However, the use of multispeed lattices 
(lattices with more then two non-zero speeds in one- 
dimension) using higher-order Gauss-Hermite quadra- 
ture is precluded because the zeros of the Hermite poly- 
nomials, for degree greater than 3, can not be aligned 
with any regular lattice. 

An alternative to Gauss-Hermite is to use Newton- 
Cotes quadrature which has the advantage that the ve- 
locities can be aligned with a lattice before determin- 
ing the quadrature weights. Multispeed lattices can be 
readily constructed. For example, suppose we want to 
employ 5 symmetric quadrature nodes {vq, . . . jV^} := 
{0,±aC,±bC}. Then, w, = CW^^ exp {v(/C^), where 



4a262 _ 



Wo 



Wi=W2^ 



2&2 



4a262 
8(52 -a2)a2' 



W3 = Wi = 



2a2 



8(a2 - 62)52 



For five distinct nodes we should have < a < 5 but 
there should be further restrictions on a and b to ensure 
all weights are positive. 

A selection of different quadrature rules are collected in 
TableUalong with their respective lattice "sound speeds" . 



IV. NUMERICAL EXPERIMENT 

The ID shock tube for a compressible isothermal fluid 
is a standard benchmark test for hydrodynamic codes. 
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FIG. 3: Density profile of a ID isothermal shock 
tube simulation after lOOv^/cs time steps using LBGK 
(solid line) and the coupled steps (|14|l (dashed line) 
with the lattice a) {vi,Wi,Cs) = ({0, ±1}, {§, |}, l/\/3); 
b) iv„W„c,) = ({0,±l,±2},{^,^,^},l/y2); c) 
{v,,W.,c,) = ({0, ±1, ±2, ±3}, i, ^, 1) 

We will fix the kinematic viscosity of the fluid a.t ly = 
10~^. Our computational domain will be the interval 
[0, 1] and we discretize this interval with 801 uniformly 
spaced lattice sites. We choose the initial density ratio as 
1:2 so that for x < 400 we set n = 1.0 else we set n = 0.5. 

In all of our simulations we use a lattice with unit 
spacing, a unit time step and we have chosen to consider 
the three-velocity Gauss-Hermite model, as well as five- 
and seven- velocity Newton-Cotes models (see Fig. [3] for 
the specific lattices details). 

The standard three-velocity LBGK discretization ([13]) 
is violently oscillatory in the immediate neighbourhood 
of the shock, the effects of which extend over much of 
the domain. The situation is worse for LBGK with five- 
and seven-velocities; here the oscillations do not remain 
bounded (we always implement the positivity rule [1, 0| 
that prohibits negative densities). The scheme quickly 
becomes meaningless. In comparison, we find that the 



proposed scheme of coupled steps p4)) gives stable be- 
haviour on all tested lattices, by which we mean no blow- 
up and spurious post-shock oscillations are observed. In 
all cases we do observe a small deviation near the shock 
front which may be an unavoidable Gibbs effect. The am- 
plitude of the Gibbs effect in our experiments decreases 
by increasing the degree of approximation. 

V. CONCLUSIONS 

By considering a realisation in which the LBM is 
wholly generated by free-flight and entropic involution, 
we have suggested a simple stabilisation procedure of cou- 
pled steps ^ which retains the second-order accuracy of 
the method without artificial viscosity at that order. 

The procedure effectively constitutes no additional 
computational cost and can be implemented in any ex- 
isting LBM code with just a few changes. 

We have demonstrated that the oscillatory pattern of 
the LBM in the vicinity of shocks is not due to a lack 
of artificial diffusivity and are not pertinent to proper 
lattice Boltzmann schemes. On the other hand, there 
exists a Gibbs effect that we cannot fully suppress with- 
out artificial viscosity or special ad hoc approximation 
schemes. 

The role of Gauss-Hermite quadratures in LBM con- 
structions is now overestimated (here we fully agree 
with [3| ) . Other quadratures require slightly more points 
for the same degree of accuracy, but can be realised on 
integer nodes (and, more generally, are flexible with re- 
gards to the choice of nodes), which is more convenient 
for LBM needs. 

The five-velocity lattice {0, ±1, ±2} can produce a sta- 
ble LBGK based method, as can many other "suspi- 
cious" lattices which were previously discussed by many 
authors as intrinsically unstable. These results can be 
immediately generalised onto high dimensional lattices 
(A'^D-lattices, cubatures and equilibria as products of ID 
lattices, quadratures and equilibria). Finally, with the 
advent of stable multispeed lattice formulations comes 
the prospect of stable nonisothermic LBM realisations. 
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